## pval_cutoff: 0.05
## lfc_cutoff: 1
## low_counts_cutoff: 10

General statistics

# Number of samples
length(counts_data)
## [1] 6
# Number of genes
nrow(counts_data)
## [1] 55487
# Total counts
colSums(counts_data)
## SRR13535276 SRR13535278 SRR13535280 SRR13535288 SRR13535290 SRR13535292 
##     3107284     2321609     3701956     7929174     6330905     3686532

Create DDS objects

# Create DESeqDataSet object
dds <- get_DESeqDataSet_obj(counts_data, ~ treatment)
## [1] TRUE
## [1] TRUE
## [1] "DESeqDataSet object of length 55487 with 0 metadata columns"
## [1] "DESeqDataSet object of length 15188 with 0 metadata columns"
colData(dds)
## DataFrame with 6 rows and 25 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider               DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name                   source_name   SRA Study
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>                    <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>                   <character> <character>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating myoblasts   SRP303354
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating myoblasts   SRP303354
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating myoblasts   SRP303354
## SRR13535288     RNA-Seq        300 12863728500 PRJNA694971 SAMN17587373 5128876770         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943372         C          GSM5043450 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043450 C2C12 proliferating myoblasts   SRP303354
## SRR13535290     RNA-Seq        300 12849825300 PRJNA694971 SAMN17587371 5136077921         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943374         C          GSM5043454 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043454 C2C12 proliferating myoblasts   SRP303354
## SRR13535292     RNA-Seq        300 10569142200 PRJNA694971 SAMN17587369 4229018065         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943376         C          GSM5043457 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043457 C2C12 proliferating myoblasts   SRP303354

Sample-to-sample comparisons

# Transform data (blinded rlog)
rld <- get_transformed_data(dds)

PCA plot

pca <- rld$pca
pca_df <- cbind(as.data.frame(colData(dds)) %>% rownames_to_column(var = 'name'), pca$x)
summary(pca)
## Importance of components:
##                            PC1     PC2     PC3     PC4    PC5      PC6
## Standard deviation     30.6369 29.5690 26.7750 23.5645 20.988 9.01e-14
## Proportion of Variance  0.2662  0.2480  0.2033  0.1575  0.125 0.00e+00
## Cumulative Proportion   0.2662  0.5142  0.7176  0.8750  1.000 1.00e+00
ggplot(pca_df, aes(x = PC1, y = PC2, color = label)) +
  geom_point() +
  geom_text(aes(label = name), position = position_nudge(y = -2), show.legend = F, size = 3) +
  scale_color_manual(values = colors_default) +
  scale_x_continuous(expand = c(0.2, 0))

Correlation heatmap

pheatmap(
  cor(rld$matrix),
  annotation_col = as.data.frame(colData(dds)) %>% select(label),
  color = brewer.pal(8, 'YlOrRd')
)

Wald test results

# DE analysis using Wald test
dds_full <- DESeq(dds)
colData(dds_full)
## DataFrame with 6 rows and 26 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider               DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name                   source_name   SRA Study        sizeFactor
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>                    <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>                   <character> <character>         <numeric>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating myoblasts   SRP303354 0.679517795445061
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating myoblasts   SRP303354 0.949141263638101
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating myoblasts   SRP303354 0.775037235318696
## SRR13535288     RNA-Seq        300 12863728500 PRJNA694971 SAMN17587373 5128876770         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943372         C          GSM5043450 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043450 C2C12 proliferating myoblasts   SRP303354  1.90346170413397
## SRR13535290     RNA-Seq        300 12849825300 PRJNA694971 SAMN17587371 5136077921         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943374         C          GSM5043454 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043454 C2C12 proliferating myoblasts   SRP303354  1.29343465761628
## SRR13535292     RNA-Seq        300 10569142200 PRJNA694971 SAMN17587369 4229018065         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943376         C          GSM5043457 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA    in space with gravity  2021-09-09  GSM5043457 C2C12 proliferating myoblasts   SRP303354 0.820640543694503
# Wald test results
res <- results(
  dds_full,
  contrast = c('treatment', condition, control),
  alpha = pval_cutoff
)
res
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 15188 rows and 6 columns
##                            baseMean     log2FoldChange             lfcSE               stat              pvalue              padj
##                           <numeric>          <numeric>         <numeric>          <numeric>           <numeric>         <numeric>
## ENSMUSG00000098104 5.46467580943419   1.19222272527557 0.980392825585785   1.21606635030526   0.223959647695148                NA
## ENSMUSG00000103922 2.19054963487128  -0.26440015750769  2.00489254415536  -0.13187747057989   0.895081208841364                NA
## ENSMUSG00000033845 170.961551920931 -0.410659544303146 0.542058807911242 -0.757592235952355   0.448695129575756 0.897546598917915
## ENSMUSG00000102275 2.42482713116482  0.445772822746135  1.59993032331677  0.278620147546185   0.780536349253174                NA
## ENSMUSG00000025903 145.031917097888 -0.144158717997538 0.305940505671955 -0.471198534763856   0.637498964963172 0.941448632245291
## ...                             ...                ...               ...                ...                 ...               ...
## ENSMUSG00000061654 89.8396453928715  -5.63949816718769  3.75874932034961  -1.50036559678397                  NA                NA
## ENSMUSG00000079834 57.8910348496526 -0.191889534789484 0.531060354122093 -0.361332818953696   0.717850662655292 0.959205276501084
## ENSMUSG00000095041 282.407505680812 -0.142150831575804 0.680850436541708 -0.208784226236002   0.834616685616294 0.980565720274868
## ENSMUSG00000063897 38.5880238283551  0.249798899406184 0.436318288783867  0.572515307809899   0.566972918209039 0.928149881015563
## ENSMUSG00000095742 8.84404223926087   2.66871532821583 0.872879538070249   3.05736955882342 0.00223288800281555  0.13808205749417
mcols(res)
## DataFrame with 6 rows and 2 columns
##                        type                               description
##                 <character>                               <character>
## baseMean       intermediate mean of normalized counts for all samples
## log2FoldChange      results  log2 fold change (MLE): treatment A vs C
## lfcSE               results          standard error: treatment A vs C
## stat                results          Wald statistic: treatment A vs C
## pvalue              results       Wald test p-value: treatment A vs C
## padj                results                      BH adjusted p-values
summary(res)
## 
## out of 15188 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 9, 0.059%
## LFC < 0 (down)     : 72, 0.47%
## outliers [1]       : 179, 1.2%
## low counts [2]     : 3527, 23%
## (mean count < 8)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotDispEsts(dds_full)

Summary details

# Upregulated genes (LFC > 0)
res_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res[which(is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 179 rows and 6 columns
##                            baseMean    log2FoldChange            lfcSE              stat    pvalue      padj
##                           <numeric>         <numeric>        <numeric>         <numeric> <numeric> <numeric>
## ENSMUSG00000103509 9.61464501255551 -4.46237357619251 2.56029953073349 -1.74291075033479        NA        NA
## ENSMUSG00000079554 38.1513589616965 -6.55179021966103 2.04362077425889 -3.20597162750854        NA        NA
## ENSMUSG00000085842 32.4579801015249   3.9244091185441  2.2620877809193  1.73486155207878        NA        NA
## ENSMUSG00000103553  13.011590868448 -4.72236490360827  2.9074692870391 -1.62421832782881        NA        NA
## ENSMUSG00000102425 25.5983581230653  6.42741446982686 2.37852118878299  2.70227337058769        NA        NA
## ...                             ...               ...              ...               ...       ...       ...
## ENSMUSG00000024867 25.1587588691625 -1.33805691056957 1.33535206509993 -1.00202556729445        NA        NA
## ENSMUSG00000117704 65.4156140694741 -4.86080355210405 1.98458148895025 -2.44928393173474        NA        NA
## ENSMUSG00000025089 56.0189687241583 -1.68541394912476 1.37461618678161 -1.22609784849895        NA        NA
## ENSMUSG00000048029 24.8195716674746 -6.69368791055847  3.9082810108819 -1.71269360926738        NA        NA
## ENSMUSG00000061654 89.8396453928715 -5.63949816718769 3.75874932034961 -1.50036559678397        NA        NA
# Low counts (only padj is NA)
res[which(is.na(res$padj) & !is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 3527 rows and 6 columns
##                            baseMean      log2FoldChange             lfcSE                stat            pvalue      padj
##                           <numeric>           <numeric>         <numeric>           <numeric>         <numeric> <numeric>
## ENSMUSG00000098104 5.46467580943419    1.19222272527557 0.980392825585785    1.21606635030526 0.223959647695148        NA
## ENSMUSG00000103922 2.19054963487128   -0.26440015750769  2.00489254415536   -0.13187747057989 0.895081208841364        NA
## ENSMUSG00000102275 2.42482713116482   0.445772822746135  1.59993032331677   0.278620147546185 0.780536349253174        NA
## ENSMUSG00000103280 3.07359630335353  -0.631912445629046   1.2048909125565  -0.524456147061706 0.599961313256629        NA
## ENSMUSG00000033740 2.71686422103861   -2.68143373680077  2.56608325087095    -1.0449519655649 0.296045170568551        NA
## ...                             ...                 ...               ...                 ...               ...       ...
## ENSMUSG00000064342 5.31912308829307  -0.322504137069612   1.2944673204709  -0.249140423994862 0.803252161252163        NA
## ENSMUSG00000064344 5.55025123689903  -0.733353272504932  1.32945753892904  -0.551618424079714 0.581209810777902        NA
## ENSMUSG00000064349 4.21295092586664 -0.0813580581645715  1.11002866545836 -0.0732936551066244 0.941572440613444        NA
## ENSMUSG00000064358 2.53259528400106    1.66589014336128  1.34176983185559    1.24156178191713 0.214398289974574        NA
## ENSMUSG00000064369  6.6742783733321    0.86939315845224 0.980303987049145   0.886860779857927 0.375153859560243        NA

Shrunken LFC results

plotMA(res)

# Shrunken LFC results
res_shrunken <- lfcShrink(
  dds_full,
  coef = str_c('treatment_', condition, '_vs_', control),
  type = 'apeglm'
)
res_shrunken
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 15188 rows and 5 columns
##                            baseMean       log2FoldChange             lfcSE              pvalue              padj
##                           <numeric>            <numeric>         <numeric>           <numeric>         <numeric>
## ENSMUSG00000098104 5.46467580943419   0.0448465000467892 0.196771185236823   0.223959647695148                NA
## ENSMUSG00000103922 2.19054963487128 -0.00240059215008312 0.191823311959347   0.895081208841364                NA
## ENSMUSG00000033845 170.961551920931  -0.0468361261370054 0.188929482787428   0.448695129575756 0.897546598917915
## ENSMUSG00000102275 2.42482713116482  0.00639198418228839 0.191481742427142   0.780536349253174                NA
## ENSMUSG00000025903 145.031917097888   -0.041440107418986 0.167089292710563   0.637498964963172 0.941448632245291
## ...                             ...                  ...               ...                 ...               ...
## ENSMUSG00000061654 89.8396453928715 -0.00740276489536808 0.192885045874907                  NA                NA
## ENSMUSG00000079834 57.8910348496526  -0.0221472788783834 0.182797420790909   0.717850662655292 0.959205276501084
## ENSMUSG00000095041 282.407505680812  -0.0102447131743499 0.185774300288715   0.834616685616294 0.980565720274868
## ENSMUSG00000063897 38.5880238283551   0.0416776527425669  0.18152022206408   0.566972918209039 0.928149881015563
## ENSMUSG00000095742 8.84404223926087      1.8996514517983  1.10322539468332 0.00223288800281555  0.13808205749417
plotMA(res_shrunken)

mcols(res_shrunken)
## DataFrame with 5 rows and 2 columns
##                        type                               description
##                 <character>                               <character>
## baseMean       intermediate mean of normalized counts for all samples
## log2FoldChange      results  log2 fold change (MAP): treatment A vs C
## lfcSE               results            posterior SD: treatment A vs C
## pvalue              results       Wald test p-value: treatment A vs C
## padj                results                      BH adjusted p-values
summary(res_shrunken, alpha = pval_cutoff)
## 
## out of 15188 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 9, 0.059%
## LFC < 0 (down)     : 72, 0.47%
## outliers [1]       : 179, 1.2%
## low counts [2]     : 3527, 23%
## (mean count < 8)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Summary details

# Upregulated genes (LFC > 0)
res_shrunken_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_shrunken_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res_shrunken[which(is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 179 rows and 5 columns
##                            baseMean       log2FoldChange             lfcSE    pvalue      padj
##                           <numeric>            <numeric>         <numeric> <numeric> <numeric>
## ENSMUSG00000103509 9.61464501255551  -0.0173337016313202 0.193722511383804        NA        NA
## ENSMUSG00000079554 38.1513589616965  -0.0340687331700345 0.197165077497607        NA        NA
## ENSMUSG00000085842 32.4579801015249   0.0191760352566012 0.193937547808487        NA        NA
## ENSMUSG00000103553  13.011590868448  -0.0129806187252493 0.193266149925685        NA        NA
## ENSMUSG00000102425 25.5983581230653   0.0226261282085071 0.194638666145371        NA        NA
## ...                             ...                  ...               ...       ...       ...
## ENSMUSG00000024867 25.1587588691625   -0.025926542768652 0.193588517717408        NA        NA
## ENSMUSG00000117704 65.4156140694741  -0.0268054351155857 0.195361376111712        NA        NA
## ENSMUSG00000025089 56.0189687241583  -0.0297509352312334  0.19468689987549        NA        NA
## ENSMUSG00000048029 24.8195716674746 -0.00752951146424761 0.192900745867621        NA        NA
## ENSMUSG00000061654 89.8396453928715 -0.00740276489536808 0.192885045874907        NA        NA
# Low counts (only padj is NA)
res_shrunken[which(is.na(res_shrunken$padj) & !is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 3527 rows and 5 columns
##                            baseMean       log2FoldChange             lfcSE            pvalue      padj
##                           <numeric>            <numeric>         <numeric>         <numeric> <numeric>
## ENSMUSG00000098104 5.46467580943419   0.0448465000467892 0.196771185236823 0.223959647695148        NA
## ENSMUSG00000103922 2.19054963487128 -0.00240059215008312 0.191823311959347 0.895081208841364        NA
## ENSMUSG00000102275 2.42482713116482  0.00639198418228839 0.191481742427142 0.780536349253174        NA
## ENSMUSG00000103280 3.07359630335353  -0.0159916365124215 0.191196494106891 0.599961313256629        NA
## ENSMUSG00000033740 2.71686422103861  -0.0134170705778806 0.193043718124969 0.296045170568551        NA
## ...                             ...                  ...               ...               ...       ...
## ENSMUSG00000064342 5.31912308829307 -0.00681405379806286 0.190793669449769 0.803252161252163        NA
## ENSMUSG00000064344 5.55025123689903   -0.014632685762997 0.191601765157723 0.581209810777902        NA
## ENSMUSG00000064349 4.21295092586664 -0.00230977701249731 0.189870502471388 0.941572440613444        NA
## ENSMUSG00000064358 2.53259528400106   0.0339351464032784  0.19522701121723 0.214398289974574        NA
## ENSMUSG00000064369  6.6742783733321    0.032470491443182 0.193120576279053 0.375153859560243        NA

Visualizing results

Heatmaps

# Plot normalized counts (z-scores)
pheatmap(counts_sig_norm[2:7], 
         color = brewer.pal(8, 'YlOrRd'), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         scale = 'row',
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts
pheatmap(counts_sig_log[2:7], 
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts (top 24 DE genes)
pheatmap((counts_sig_log %>% filter(ensembl_gene_id %in% res_sig_df$ensembl_gene_id[1:24]))[2:7],
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label), 
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

Volcano plots

# Unshrunken LFC
res_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

# Shrunken LFC
res_shrunken_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

GSEA (all)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

GSEA (DE)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

System Info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] fgsea_1.12.0                Rcpp_1.0.3                  RColorBrewer_1.1-2          pheatmap_1.0.12             DESeq2_1.26.0               SummarizedExperiment_1.16.1 DelayedArray_0.12.3         BiocParallel_1.20.1         matrixStats_0.57.0          Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.1         IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0         scales_1.1.1                forcats_0.4.0               stringr_1.4.0               dplyr_1.0.2                 purrr_0.3.3                 readr_1.3.1                 tidyr_1.0.0                 tibble_3.1.0                ggplot2_3.3.3               tidyverse_1.2.1            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       ellipsis_0.3.0         htmlTable_1.13.3       XVector_0.26.0         base64enc_0.1-3        rstudioapi_0.10        farver_2.1.0           bit64_0.9-7            mvtnorm_1.1-1          apeglm_1.8.0           AnnotationDbi_1.48.0   fansi_0.4.0            lubridate_1.7.4        xml2_1.2.2             splines_3.6.3          geneplotter_1.64.0     knitr_1.25             Formula_1.2-3          jsonlite_1.6           broom_0.7.5            annotate_1.64.0        cluster_2.1.0          png_0.1-7              compiler_3.6.3         httr_1.4.1             backports_1.1.5        assertthat_0.2.1       Matrix_1.2-18          cli_1.1.0              acepack_1.4.1          htmltools_0.5.1.1      tools_3.6.3            coda_0.19-3            gtable_0.3.0           glue_1.4.2             GenomeInfoDbData_1.2.2 fastmatch_1.1-0        bbmle_1.0.23.1         cellranger_1.1.0       jquerylib_0.1.3        vctrs_0.3.4            xfun_0.22              rvest_0.3.5            lifecycle_0.2.0        XML_3.99-0.3           MASS_7.3-51.5          zlibbioc_1.32.0        hms_0.5.2              yaml_2.2.0             memoise_1.1.0          gridExtra_2.3          emdbook_1.3.12         sass_0.3.1             bdsmatrix_1.3-4        rpart_4.1-15           latticeExtra_0.6-29    stringi_1.4.3          RSQLite_2.2.1          genefilter_1.68.0      checkmate_1.9.4        rlang_0.4.8            pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14          lattice_0.20-38        labeling_0.3           htmlwidgets_1.5.1      bit_1.1-15.1           tidyselect_1.1.0       plyr_1.8.4             magrittr_1.5           R6_2.4.0               generics_0.0.2         Hmisc_4.3-0            DBI_1.1.0              pillar_1.5.1           haven_2.2.0            foreign_0.8-75         withr_2.1.2            survival_3.1-8         RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.5           crayon_1.3.4           utf8_1.1.4             rmarkdown_2.7          jpeg_0.1-8.1           locfit_1.5-9.4         grid_3.6.3             readxl_1.3.1           data.table_1.13.6      blob_1.2.1             digest_0.6.27          xtable_1.8-4           numDeriv_2016.8-1.1    munsell_0.5.0          bslib_0.2.4